{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Pitchup at $k = 0.7$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using PotentialFlow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Plots\n",
    "\n",
    "clibrary(:colorbrewer)\n",
    "default(clim = (-0.05, 0.05), markerstrokealpha = 0, markersize = 3, grid = false, \n",
    "     legend = false, colorbar = :right, colorbar_title = \"\\$\\\\Gamma\\$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the beginning of every time-step, we first determine compute the plate kinematics, then the bound vortex sheet strength required to satisfy the no-flow-through condition, and finally velocity of all vortex elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_ẋ! (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_ẋ!(ẋ, x, t)\n",
    "    plate, ambient_sys = x\n",
    "    motion = ẋ[1]\n",
    "    motion.ċ, motion.c̈, motion.α̇ = motion.kin(t)\n",
    "    \n",
    "    Plates.enforce_no_flow_through!(plate, motion, ambient_sys, t)\n",
    "    reset_velocity!(ẋ, x)\n",
    "    self_induce_velocity!(ẋ, x, t)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Once we have advected all the vortex elements, we release new segments of the vortex sheet from the edges of the plate.\n",
    "Since we are dealing with a low angle of attack case, only the trailing edge will be forced to satisfy the Kutta condition.\n",
    "Vorticity released from the leading edge will instead follow the leading edge suction parameter (LESP) criteria described in [Ramesh et al. 2014](https://doi.org/10.1017/jfm.2014.297)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "shed_new_vorticity! (generic function with 3 methods)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function shed_new_vorticity!(blobs, plate, motion, t, lesp = 0.0, tesp = 0.0)\n",
    "    z₊ = (blobs[end-1].z + 2plate.zs[end])/3\n",
    "    z₋ = (blobs[end].z + 2plate.zs[1])/3\n",
    "    \n",
    "    blob₊ = Vortex.Blob(z₊, 1.0, δ)\n",
    "    blob₋ = Vortex.Blob(z₋, 1.0, δ)\n",
    "    Plates.enforce_no_flow_through!(plate, motion, blobs, t)\n",
    "    \n",
    "    Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, blob₊, blob₋, t, lesp, tesp);\n",
    "    \n",
    "    push!(blobs, Vortex.Blob(z₊, Γ₊, blobs[1].δ), Vortex.Blob(z₋, Γ₋, blobs[1].δ))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Motion Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We start with a plate at around 5ᵒ that smoothly pitches up to 90ᵒ, while moving forward at a constant velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "L = 1.0 # length of the plate\n",
    "\n",
    "pitchup = let U₀ = 1,    # Translation velocity\n",
    "              a  = 0.5,  # Pitching axis (leading edge)\n",
    "              K = 0.7,   # Pitch rate\n",
    "              α₀ = 0.0,  # Nominal starting angle\n",
    "              t₀ = 0.0,  # Nominal starting time\n",
    "              Δα = π/2,  # Total angle change\n",
    "              ramp = Plates.RigidBodyMotions.EldredgeRamp(6) # Smooth ramp function\n",
    "    \n",
    "    Plates.RigidBodyMotions.Pitchup(U₀, a, K, α₀, t₀, Δα, ramp)\n",
    "end;\n",
    "\n",
    "motion = Plates.RigidBodyMotion(pitchup);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Discretization and Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "N = 128\n",
    "plate = Plate(N, L, zero(Complex128), pitchup.α(0))\n",
    "\n",
    "Δt = 1e-2\n",
    "\n",
    "δ = 0.01\n",
    "lesp = 0.2\n",
    "tesp = 0.0\n",
    "\n",
    "Δz₀ = im*3Δt*exp(im*plate.α)\n",
    "z₋, z₊ = plate.zs[[1,N]]\n",
    "\n",
    "blobs = Vortex.Blob.(Δz₀ .+ [z₊, z₋], 1.0, δ)\n",
    "\n",
    "Plates.enforce_no_flow_through!(plate, motion, (), 0)\n",
    "Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, blobs[1], blobs[2], 0.0, lesp, tesp);\n",
    "\n",
    "blobs = Vortex.Blob.(Δz₀ .+ [z₊, z₋], [Γ₊, Γ₋], δ)\n",
    "\n",
    "sys₀ = (plate, blobs)\n",
    "\n",
    "sys = deepcopy(sys₀)\n",
    "sys₊ = deepcopy(sys₀) # Used for storage during time-marching\n",
    "ẋs = (motion, allocate_velocity(blobs))\n",
    "\n",
    "forces = Complex128[];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Marching"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "T = 0.0:Δt:8.0\n",
    "\n",
    "for t in T\n",
    "    \n",
    "    plate, ambient_ω = sys\n",
    "    motion, ambient_u = ẋs\n",
    "\n",
    "    resize!(sys₊[2], length(sys[2]))\n",
    "    forward_euler!(sys₊, sys, t, Δt, compute_ẋ!, advect!, ẋs)\n",
    "\n",
    "    # The force requires information about the motion of the plate,\n",
    "    # the strength, location, and velocity of the ambient vortex elements,\n",
    "    # as well as the vorticity flux from the plate edges\n",
    "    push!(forces, Plates.force(plate, motion, ambient_ω, ambient_u,\n",
    "            (ambient_ω[end-1], ambient_ω[end]), Δt))\n",
    "\n",
    "    sys, sys₊ = sys₊, sys\n",
    "    shed_new_vorticity!(sys[2], sys[1], ẋs[1], t, lesp, tesp)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sys, color = :RdBu_r, ratio = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(T, [-2real.(forces) 2imag.(forces)], layout = (2,1), linewidth = 2,\n",
    "     xlabel = \"Convective Time\", ylabel = [\"\\$C_D\\$\" \"\\$C_L\\$\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Alternate ramp functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The `Motion` type has four fields\n",
    "- `ċ` is the current centroid velocity\n",
    "- `c̈` is the current centroid acceleration\n",
    "- `α̇` is the current angular velocity\n",
    "- `kin` is a function that returns a tuple `(ċ, c̈, α̇)` given an input time\n",
    "\n",
    "Technically the first three fields are redundunt with the `kin` function, but we don't want to have to carrying around the current time everywhere in our code, so we just use the first three fields as cached values.\n",
    "In this particular test case, the `kin` function is a `Pitchup` object.  `Pitchup` takes in the tranlational velocity, pitch rate, etc., and also a `ramp` field.  It expects the `ramp` field to be a function that takes in time and returns a smoothed ramp function (with unit nominal slope).  It then composes `ramp` with a time-shifted and negated version of `ramp` to create the pitching profile.\n",
    "\n",
    "The previous example used a `logcosh` ramp function.  The following example will use a power series based ramp function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "pitchup = let U₀ = 1,    # Translation velocity\n",
    "              a  = 0.5,  # Pitching axis (leading edge)\n",
    "              K = 0.7,   # Pitch rate\n",
    "              α₀ = 0.0,  # Nominal starting angle\n",
    "              t₀ = 0.0,  # Nominal starting time\n",
    "              Δα = π/2,  # Total angle change\n",
    "              ramp = Plates.RigidBodyMotions.ColoniusRamp(4) # Smooth ramp function\n",
    "    \n",
    "    Plates.RigidBodyMotions.Pitchup(U₀, a, K, α₀, t₀, Δα, ramp)\n",
    "end;\n",
    "\n",
    "motion = Plates.RigidBodyMotion(pitchup);\n",
    "plate = Plate(N, L, zero(Complex128), pitchup.α(0))\n",
    "\n",
    "sys = (plate, sys₀[2])\n",
    "sys₊ = deepcopy(sys)\n",
    "forces = Complex128[];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "for t in T\n",
    "    \n",
    "    plate, ambient_ω = sys\n",
    "    motion, ambient_u = ẋs\n",
    "\n",
    "    resize!(sys₊[2], length(sys[2]))\n",
    "    forward_euler!(sys₊, sys, t, Δt, compute_ẋ!, advect!, ẋs)\n",
    "\n",
    "    # The force requires information about the motion of the plate,\n",
    "    # the strength, location, and velocity of the ambient vortex elements,\n",
    "    # as well as the vorticity flux from the plate edges\n",
    "    push!(forces, Plates.force(plate, motion, ambient_ω, ambient_u,\n",
    "            (ambient_ω[end-1], ambient_ω[end]), Δt))\n",
    "\n",
    "    sys, sys₊ = sys₊, sys\n",
    "    shed_new_vorticity!(sys[2], sys[1], ẋs[1], t, lesp, tesp)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sys, color = :RdBu_r, ratio = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(T, [-2real.(forces) 2imag.(forces)], layout = (2,1), linewidth = 2,\n",
    "     xlabel = \"Convective Time\", ylabel = [\"\\$C_D\\$\" \"\\$C_L\\$\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
